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Abstract 


Three random number generators, which produce Gaussian white noise sequences, were 
compared to assess their suitability in aircraft dynamic modeling applications. The first gen- 
erator considered was the MATLAB® implementation of the Mersenne-Twister algorithm. 
The second generator was a website called Random.org, which processes atmospheric noise 
measured using radios to create the random numbers. The third generator was based on 
synthesis of the Fourier series, where the random number sequences are constructed from 
prescribed amplitude and phase spectra. A total of 200 sequences, each having 601 random 
numbers, for each generator were collected and analyzed in terms of the mean, variance, nor- 
mality, autocorrelation, and power spectral density. These sequences were then applied to 
two problems in aircraft dynamic modeling, namely estimating stability and control deriva- 
tives from simulated onboard sensor data, and simulating flight in atmospheric turbulence. 
In general, each random number generator had good performance and is well-suited for air- 
craft dynamic modeling applications. Specific strengths and weaknesses of each generator 
are discussed. For Monte Carlo simulation, the Fourier synthesis method is recommended 
because it most accurately and consistently approximated Gaussian white noise and can be 
implemented with reasonable computational effort. 
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Nomenclature 


vertical accelerometer output, g 

Fourier series amplitude 

covariance and variance 

expectation operator 

frequency, Hz 

Nyquist frequency, = 1/2At Hz 
one-sided power spectral density 
gravitational acceleration, = 32.174 ft/s? 
sample index variable 

cost function 

imaginary number, = /—1 

frequency index variable 

lag index variable 

number of frequencies 

pitching moment stability and control derivatives 
number of samples 

normal distribution having zero mean and unit variance 
Probability 

pitch rate, rad/s 

coefficient of determination 
autocorrelation function 

Laplace operator 

standard error 

record length, s 

time, s 

airspeed, ft/s 

noise sequence 

vertical body-fixed velocity, ft/s 
regressor matrix 

vertical body-axis force stability and control derivatives 
dependent modeling variable 


angle of attack, rad 
perturbation quantity 

sampling period, s 

elevator deflection, deg 

model parameters 

pairwise correlation coefficient 
standard deviation 

cumulative distribution function 
Fourier series phase angle, rad 
angular frequency, rad/s 


Superscripts 


Subscripts 


complex conjugate 
time derivative 
estimated value 
mean value 
inverse 


trim value 
gust 


1 Introduction 


Random number sequences are often used in aircraft dynamic modeling applications 
to approximate Gaussian, band-limited, white noise sequences having the following ideal 
characteristics: 


e Zero mean e Uniform power spectrum 
e Unit variance e No serial correlations 
e Gaussian/normal distribution e Finite energy 


For example, this type of noise is often added to simulated output signals to mimic 
sensor measurement noise. White noise can also be used to drive a dynamic system to 
simulate sources of process noise. Shaping or coloring filters use white noise as input to 
produce random signals with a specified power spectral density, for example to simulate 
neglected deterministic dynamics [1,2] or stochastic turbulence [3]. Drift in computer clocks, 
such as in GPS receivers, or drift in sensor calibration parameters, for example due to 
temperature sensitivities, can be modeled as a random walk using white noise [4]. In these 
modeling applications, many different white noise sequences are used in repeated Monte 
Carlo experiments to gain confidence in the predictive capability of a model, to evaluate 
handling qualities or control law performance, or to develop and test new methods of analysis 
before application to flight test data. In identification applications, the performance of 
parameter estimation algorithms using the extended Kalman filter improve when artificial 
white noise is used in the analysis to allow the parameter estimates to vary [5]. 


Given the ubiquity of random numbers in aircraft dynamic modeling problems, it is 
important to have fast access to many random number sequences of high quality and po- 
tentially long duration. Algorithms for generating random number sequences have been 
improving since the advent of the modern computer [6]. During the 1980’s and 1990s, 
as Monte Carlo simulations of aerospace vehicles became more widespread, many methods 
for generating random numbers were developed and evaluated. More recently, the focus 
has been on generating random numbers from measurements of seemingly random physical 
processes, especially for cryptography applications. 


Despite this steady growth of research, no recent studies on evaluating random num- 
ber generators for aircraft dynamic modeling applications were found in the literature. A 
casual polling of colleagues indicated that most users had not investigated the quality of 
their random number generators, but generally felt comfortable with the results for their 
applications. On this topic of using random number generators to simulate noise, Hamming 
writes [7] 


Experience seems to indicate that the careless use of a random-number generator 
will trap the unwary many more times than at first seems reasonable, but also 
that with care random-number generators give very satisfactory results. It is a 
field for the beginner to be wary and suspicious of and he should make a few 
careful checks that all is going as he thinks it should before he plunges into the 
details of using the simulation as a working tool of design optimization. 


This report documents such an investigation. The purpose is to examine a few interesting 
random number generators, using relevant statistical tools, and report on the suitability of 
the generators for use in dynamic modeling applications. 


Three sources of random number sequences were considered. The first was the Mersenne- 
Twister algorithm implemented as the default random number generator in the MATLAB® 
software suite. The second was the website Random.org, which streams random numbers 
processed from atmospheric noise measured using radios. The third was a method involving 
the construction of random sequences through Fourier synthesis. 


The following section describes the random number generators. Results are then shown 
for one sample record, to illustrate the characteristics of these sequences using statistical 
tools commonly used for dynamic modeling. To examine the asymptotic nature of the 
generators, results using ensembles of sequences are then presented. These same random 
sequences are then applied to two example aircraft modeling applications in a Monte Carlo 
simulation framework. The first example is estimating stability and control derivatives from 
simulated onboard sensor data using linear regression. The second example is adding at- 
mospheric turbulence to an aircraft flight dynamics simulation. Concluding remarks then 
complete the report. In general, each random number generator performed well and was 
considered suitable for Monte Carlo simulation. The method of Fourier synthesis performed 
better than the other two methods in that it more consistently and more precisely approxi- 
mated Gaussian white noise. 


Software for the multisine input design and much of the analysis are available in a 
MATLAB toolbox called System IDentification Programs for AirCraft, or SIDPAC [8], which 
is associated with Ref. 9. 


2 Random Number Generators 


MATLAB / Mersenne-Twister 


The MATLAB [10] software suite has a built-in random number generator. The current 
software version (R2016b) invokes the Mersenne-Twister method [11], which was developed 
specifically for Monte Carlo simulation. The software creates a sequence of 32-bit random 
integers, which are then transformed to the open interval (0,1). The function randn.m takes 
this uniformly-distributed sequence and maps it to a zero-mean, unit-variance, normally- 
distributed sequence using the Ziggurat algorithm [12]. Other non-default settings for the 
MATLAB random number generator were not investigated because they implement legacy 
algorithms or other distributions. 


This generator was considered in this report for several reasons. Because MATLAB 
is widely used in the analysis of aerospace vehicles, its native random number generator 
is used often. Mersenne-Twister is a deterministic algorithm known as a pseudo-random 
number generator, and previously-used sequences can be regenerated or new sequences can 
be generated. This is advantageous for debugging software, but it also means that the 
sequences are not truly random, and that there is a finite period before the sequence repeats. 
However, this method can create many long sequences of random numbers quickly, which 
facilitates Monte Carlo simulation. 


Random.org 


The website Random.org [13] streams random numbers created from measurements of 
atmospheric noise using radio receivers. The radios are set to unused frequencies to register 
static noise, which is digitized as a sequence of 8-bit numbers. All bits in these numbers are 
discarded, except the least significant bit, which has a high level of entropy. Retained bits 
from many different recorded numbers are then combined to form a uniformly-distributed 
number. For a Gaussian sequence, the Box-Muller method [14] is used to transform the 
distribution. 


This random number generator was considered because it is regarded as a true random 
number generator, which should be free of the problems associated with pseudo-random 
number generators. In addition, the website contains documentation by the Trinity College, 
Dublin evaluating the quality of the results from the website. Therefore this generator 
was expected to be of high quality and useful for comparison. Furthermore, access to the 
website was free of charge, and the user interface provided enough control to facilitate a fair 
comparison with other methods, e.g., distribution type, number of significant digits, etc. 


A relatively large amount of time was needed to query the website for random num- 
bers, and then to store and recall the results for analysis. As Monte Carlo simulations 
involve repeated calculations and require relatively large amounts of time, any additional 
delay due to obtaining or handling random numbers is not desirable. Another drawback 
is that some true random number generators are susceptible to environmental factors that 
degrade performance, such as the aging of electrical components that may bias the numbers. 
Random.org guards against these factors by using multiple radios in different locations and 
monitoring statistical metrics in real time. 


Fourier Synthesis 


Lanczos and Bellai [15] expound on an earlier observation by Lanczos [16] that an ideal 
Gaussian white noise sequence is characterized by a constant amplitude spectrum and a 
random phase spectrum, and that such a sequence can be synthesized using the Fourier 
series. Specifically, given N equidistant sample times 
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and M harmonic frequencies up to but not including the Nyquist frequency 


k E 
faa fot RSL Ahsan (2) 


a Gaussian random noise sequence can be constructed using the finite Fourier series as 


v(t;) = s Cy sin (Fe re én) (3) 


k=1 


The amplitudes c, are set equal to a constant to give a flat amplitude spectrum, charac- 
teristic of white noise. Choosing this constant as ,/2/M yields a sequence with unit variance, 
as shown in Appendix A. The phase angles ¢; are drawn from a uniform distribution over 
the interval (0,277). Any method of obtaining uniformly-distributed random numbers can 
be used for this step. In this report, the MATLAB function rand.m was used to implement 
the Mersenne-Twister algorithm described in Section 2. The resulting sequences have, to 
numerically accuracy, zero mean because the bias term in the Fourier series was omitted and 
because the harmonic sinusoids had complete cycles over the time duration [0,7]. When 
the sequence is long enough that many harmonic sinusoids are present in the sequence, the 
central limit theorem suggests the sequence attains a Gaussian distribution [15]. 


This generator was considered because of its excellent potential qualities. It also retains 
the favorable properties of the pseudo-random number generator in that sequences can be 
constructed quickly and can be regenerated for software debugging or randomized for Monte 
Carlo simulation. Another benefit is that the Fourier series is used to both construct the 
sequence and to analyze its spectral character. As discussed later, this results in an excellent 
approximation of Gaussian white noise in the frequency domain. If this method is to be 
used in real time, for example with flight simulators or parameter estimation using extended 
Kalman filters, a large value of T can be selected to synthesize a noise sequence. If that 
value of time is reached, a new sequence can be constructed using another random drawing 
of phase angles. 


3 Single Sample Records 


In this section, a single sample record from each random number generator is studied 
to illustrate typical statistical characteristics of interest in dynamic modeling. The next 
section examines the asymptotic qualities of the generators. 


Standardized tests for random number generators are typically intended for uniform 
distributions or binary number sequences, and evaluate a variety of characteristics such as 
spectral content, pattern reoccurrence, total draws of 0’s and 1’s, and furthest drift incurred 
from random walk processes [17,18]. These tests are not all directly applicable or relevant 
to the random number sequences studied here. Instead, this report uses analysis tools 
commonplace in aircraft dynamic modeling applications to judge the precision in which these 
sequences approximate Gaussian white noise. In this context, a successful generator will 
closely match all the prescribed characteristics of Gaussian white noise mentioned previously. 


Each sequence contains 601 random numbers in double-precision format. The length of 
the sequence was selected to reflect typical aircraft modeling data records. Specifically, these 
sequences correspond to 12 s of data recorded at 50 Hz, which could be a maneuver involving 
the short period or dutch roll motions of an aircraft and standard onboard instrumentation. 
Each set of random numbers was generated as a Gaussian, zero-mean, unit-variance, white 
noise sequence. Figure 1 shows each of the three sample records. The plot markers are used 
consistently throughout this paper to distinguish the random number generators. Summary 
statistics for these sequences are listed in Table 1. 


Mean and Variance 


The random number sequences were intended to have a mean of zero and a variance 
of unity. For the sample records shown, the MATLAB sequence had the largest error in 
the mean, and the Random.org sequence had the largest error in the variance. The Fourier 
synthesis sequence had almost no error in these regards, by design. The small amount 
of error in the mean was due to numerical rounding errors, which could be removed by 
enforcing boundary conditions. 


Although the variances of the sample records are approximately the same, Fig. 1 shows 
that there is a greater dispersion of the data in the sequences from Random.org and Fourier 
synthesis than from MATLAB. These points are not statistical outliers, as discussed in the 
following section. The comments given here applied in general, not just for these particular 
sample records. 


Normality 


Normal probability plots show data and assumed theoretical distributions, drawn as 
straight lines [19]. These plots are useful for detecting systematic and large deviations 
from a theoretical distribution, and do not require any adjustments by the analyst, such as 
selecting bin sizes and locations for histograms. 


For a normal distribution, the transformed data are [9] 


o~ [PC] (4) 


where 


POS ( 2 5) IN (5a) 


H(z) = Gal. oat (5b) 


Here P(i) is the empirical quantile quantifying the fraction of data less than its argument, 
and ®(z) is the cumulative distribution function. When the transformed data are plotted 
against the ordered data, the slope quantifies the variance and the intercept quantifies the 
mean. The degree to which the transformed data diverge from the straight line describes 
the inaccuracy of the assumed distribution. 


Normal probability plots are shown in Fig. 2 for the sample records. Also shown with 
the data is a straight line, corresponding to the normal distribution with zero mean and 
unit variance. Agreement of these plots indicates that each sequence was approximately 
Gaussian. Coefficients of determination 


ol vu — Nv? 


2 
as viv — Nv? (6) 


listed in Table 1 indicate that fits to the theoretical distribution all had less than 0.4% 
error. There were slight deviations from the theoretical distribution in the larger quartiles 
of the data, but this is expected for distributions with small tails, such as the normal 
distribution [19]. 


As remarked earlier from viewing the sample records, it is interesting that the sequences 
from Random.org and Fourier Synthesis have a larger dispersion of data, although a similar 
variance. These data correspond to the tails of the normal probability distribution function, 
as illustrated in Fig. 2 by data in the fourth quartiles. Because these data still lie close 
to the straight line, they are consistent with a normal distribution and are not statistical 
outliers. 


Autocorrelation 


The autocorrelation of a sequence can be used to show serial correlations in the data 
and evaluate whiteness. An estimate of the autocorrelation is [9] 
1 Xa! 


t=1 


Py (L) v(iv(it)), for’ =0,1,2,...,N-1 (7) 


where / is a lag index. This estimate is biased, but is accurate for data with large N, 

which is common in aircraft modeling problems [9]. The standard error of this estimate is 
approximated as [9, 20] 

Puy (0) 
VN 


For white noise sequences, plotting the autocorrelation for various lag index values re- 
sults in a peak at | = 0, equal to the variance of the sequence, and zero at all other lag 
indices. Systematic and large deviations from this description indicate the presence of serial 
correlations in the data. In that case, the data is representative of colored noise, not white 
noise. 


s[Fuo(2)] = for 1 #0 (8) 
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Figure 3 shows the autocorrelations of the sample records, as well as two standard 
error intervals, where s[7,,(1)] = 1/601. Each sequence appears approximately white. 
Table 1 indicates that the variance estimate attained from #,,(0) was less accurate for the 
Random.org sequence than for the other two. 


There are appreciable differences between the autocorrelations of the sequences for non- 
zero lag indices. The MATLAB sequence had 0.5% of the lagged autocorrelation values 
outside the error bound, and the Random.org sequence had 2.16% of its values outside the 
error bound. While these values are in statistical agreement with the 95% confidence in- 
terval, these points could be confused with colored noise sequences or unmodeled dynamics 
in dynamic modeling problems. The Fourier synthesis sequence had all of its lagged auto- 
correlation values remain well within the error bounds. The reason for this is discussed in 
Appendix B 


Power Spectral Density 


White noise has a uniform power spectral density (PSD), which can be computed as 


. Di dss: . 

Gov (Jun) = Fv" (Jwn)vGwe) (9) 
where wy, = 2nf,, v(jw,) is computed from the finite Fourier transform of the sequence 
(9, 21], and v*(jw,) denotes the complex conjugate of v(jw,). The transform was evalu- 
ated using the chirp z-transform at the frequencies given by Eq. (2), which were the same 
frequencies used in the method of Fourier synthesis to create the noise sequence. 


Figure 4 shows the PSD for each of the three sample records. These plots are multiplied 
by the Nyquist frequency so that the values shown approximate the variance of the signal. 
In addition, the means of these data (averaged over frequency) are shown, as well as two 
standard deviation intervals about those mean values. On average, each of the three sample 
records approximated unit variance white noise well in this regard; however, the performance 
of the Fourier synthesis method was clearly superior. This was a result of designing the 
random sequences in the frequency domain to have a uniform amplitude spectrum. 


The variance of the amplitude spectrum was significantly lower for the sample record 
generated using Fourier synthesis. This was in part due to its frequency-based design, and 
also because the same harmonic frequencies used to construct the sequence were used in 
the Fourier transform to analyze it. If, for example, the frequency resolution of the chirp 
z-transform were increased by a factor of 2, more scatter would be seen in the PSD of the 
Fourier synthesis data due to spectral leakage caused by using a finite duration of data. 
This was investigated, and the mean over frequency decreased to 0.9968 and the standard 
deviation increased to 0.5412. Comparison with the other methods in Table 1 shows this is 
still a better approximation of white noise than the sample records obtained from MATLAB 
and Random.org, from this perspective. One may attempt to improve the quality of the 
noise sequence by using a finer spacing of frequencies in its construction, i.e., a frequency 
resolution finer than 1/T in Eq. (2), but this forfeits the orthogonality property of the 
sinusoid series and introduces serial correlations into the sequence. Instead, a longer noise 
sequence could be designed and truncated to give the same result using a finer frequency 
resolution. 
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4 Ensembles of Sample Records 


The results in the previous section examined a single sample record from each random 
number generator. Analysis of ensembles of sequences are discussed in this section to ex- 
amine the asymptotic nature of the random number generators. In total, 200 different 
sequences of 601 random numbers were collected from each random number generator. The 
MATLAB Mersenne-Twister sequences and the Fourier synthesis sequences were generated 
at the time of analysis. The sequences from Random.org were queried thirty sequences at 
a time over seven days, and at different times of day, to remain within the allowed website 
quota and to attempt to remove any correlations due to time of day, weather, etc. 


The results of the analysis generally followed the same trends as those exhibited from 
the single sample records discussed in the previous section. Summary statistics are listed in 
Table 2. The words “average” and “scatter” are used here to describe the mean and standard 
deviation of the Monte Carlo ensembles, to distinguish them from the words “mean” and 
“variance,” which were reserved for the evaluation of the individual sequences. 


In addition, the pairwise correlation coefficient between two sequences v1 (t;) and v2(t;) 


N 
De [v1 (ti) — D1] [v2 (ti) — 02] 
p= i=1 (10) 


N N 
> [vr (ti) — oi)” [v2(t;) — 82]? 


i=l i=1 


was investigated. This metric ranges between +1 and is a measure of similarity between two 
sequences, which could arise from periodicities in the data or similar starting seed values 
for pseudo-random number generators, or from deterministic fluctuations manifesting in the 
true random number generators. Metrics listed in Table 2 indicated this was not the case 
and the sequences were not appreciably correlated with each other. Figure 5 shows two 
sequences plotted against each other for each generator. Because these data do not form a 
straight line, they visually indicate a low correlation. Each block of thirty sequences from 
Random.org showed similar values of correlation as the full ensemble. 


The MATLAB and Random.org generators achieve a normal distribution through an 
explicit mapping of the random numbers. In contrast, the Fourier synthesis method only 
suggests a Gaussian distribution via the central limit theorem. To investigate the asymptotic 
convergence of the sample distribution, a separate Monte Carlo analysis was performed using 
only the Fourier synthesis method. Sequences with N = 1,2,...,601 samples were generated 
200 different times each, and the R? values for the normal probability plots were evaluated. 
Figure 6 shows the mean and two standard deviations of scatter in the results. A normal 
sample distribution was achieved relatively quickly using Fourier synthesis. On average, 
approximately 1% error was attained after using N = 100 (50 sinusoids), which represents 
2s of data sampled at 50 Hz. 


Overall, all three of the random number generators performed well in approximating 
Gaussian white noise. The method of Fourier synthesis most consistently simulated the 
desired characteristics per the statistics shown. The Fourier synthesis method also efficiently 
attains a Gaussian distribution. 
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5 Application to Parameter Estimation 


One use of measured flight test data is to estimate aerodynamic stability and control 
derivatives. These estimates can be used in a variety of applications, such as providing phys- 
ical insight into the aircraft, constructing flight dynamics simulations, developing feedback 
control laws, evaluating handling qualities, or computing performance metrics. 


In the development and testing of a parameter estimation method, a simulation model 
is often used first before applying the estimator to measured flight test data. As confidence 
in the analysis is gained, various sources of error can be included to assess performance and 
robustness in realistic conditions. It is therefore important to accurately replicate the effects 
expected in measured flight test data. 


To investigate the effect of different random number generators on parameter estimation, 
a simple example was investigated. The simulation model was the short-period dynamics of 
an aircraft about a reference flight condition, given in state-space form as [9,22] 


&])_ [| Z 14+2Z,]f[ Aa Zs, 

[a l= [ae “he | [se + | a, | 44 aa 
Aa 1 0 x 0 
Aq | = 0 i Pa 0 Nbc (11b) 
Aa, az, eZ, q Ma Zs, 


where a is the angle of attack, q is the pitch rate, 6, is the elevator deflection, and a, is 
the vertical acceleration at the center of mass. The goal here is to estimate values of the 
stability and control derivatives, such as Z, and M;5,, from measurements of the input and 
output data. 


The model was simulated for 12 s and sampled at 50 Hz, which corresponded to N = 601 
samples, as before. The elevator input included a 10 s multisine excitation [9] with frequency 
content between 0.3 Hz and 2.1 Hz!. The model parameters for this simulation are listed 
in Table 3, and pertain to the T-2 subscale aircraft, shown in Fig. 7, in a typical flight 
condition. The frequency of the short period mode was 1.04 Hz. The time histories resulting 
from simulating Eqs. (11) are shown in Fig. 8. 


This estimation problem can be posed as the least squares problem 
z=X0+Vv (12) 


where z are the dependent variables, X are the explanatory variables, 8 are the unknown 
model parameters, and v is a Gaussian white noise sequence. Solving the least-squares cost 
function 


J(0) = ; (z — X60)" (z — X6) (13) 


yields the optimal estimate 
6 = (X7X) 'X7z (14) 


The uncertainty in this estimate is quantified by the covariance matrix 


cov(6) = (XTX) ' XT E[vv7]X (KTX)* (15) 


1Note that these multisines are also described by Eq. (3), and therefore can also be thought of as a 
truncated Fourier series expansion. 


13 


where E[vv“] is the autocorrelation matrix, formed using Eq. (7). The parameter covariance 
matrix is therefore sensitive to serial correlations in the modeling residual sequence. 


The estimation was performed twice for each Monte Carlo simulation. In the first case, 
z were the measurements Aa, and @ included the Z derivatives. In the second case, z were 
the measurements g, obtained by smooth numerical differentiation [9] of the measured gq, 
and @ included the M derivatives. In each case the matrix X consisted of Aa, Ag, and Ad, 
in column vectors. 


A Monte Carlo simulation was performed using the same 200 random noise sequences 
from each of the three random number generators. In each run, these noise sequences were 
scaled to 20% of the root mean square variation of the corresponding signal, and added to 
z. Because estimation of the Z derivatives is performed separately from the M derivatives, 
the same noise sequences can be used without correlating the estimation results. Noise was 
not added to X because these data should be smoothed before estimation to remove the 
noise. 


Model fits to one sample record of simulated Aa, and q data for each random number 
generator are shown in Figs. 9 and 10. The fits were generally good and had small residuals. 
The summary statistics for the Monte Carlo simulations are listed in Table 4. Parameter 
estimates were generally similar in value for each random number generator. Each set of 
parameter estimates were in statistical agreement with the true values, and had scatter of 
the parameters within the standard errors. The results using noise sequences generated 
using Fourier synthesis had consistently lower scatter in the estimated standard errors than 
the results using MATLAB and Random.org. The condensed scatter also indicates that the 
method of Fourier synthesis produces sequences with lower serial correlations, as evident 
from Table 2 and previous work [1, 2]. 
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6 Application to Simulation in Turbulence 


Simplified models of atmospheric turbulence are useful engineering approximations for 
evaluating handling qualities of piloted aircraft, assessing robustness of control laws and 
identification algorithms, and other applications. One such model is the continuous stochas- 
tic gust [3]. The gust history is modeled as a random variable with a specified PSD which 
is dependent upon variables including turbulence level, altitude, and airspeed. The spectral 
character of the random sequences is achieved by passing unit variance white noise through 
a coloring filter. 


In particular, the Dryden turbulence model for the body-fixed vertical velocity w has 
the filter transfer function 
Aw,(s) 2D, (1+ 2V3Lys/Vo) 


v(s) are ™Vo (14+2Lys/Vo)" oy) 


where L,, is the length scale and oy is the turbulence intensity. Parameters used for the 
simulation are listed in Table 3. Each of the white noise sequences generated was passed 
through the coloring filter to produce a vertical gust time history. For the first sample 
records, the gust sequences are shown in Fig. 11 and the PSD estimates, again multiplied by 
the Nyquist frequency, are shown in Fig. 12. The sequences from MATLAB and Random.org 
had many data points with less power than prescribed. At the higher frequencies, the 
PSD from Random.org appears to begin to diverge from the specified PSD. The sequence 
generated using Fourier synthesis matched the specified PSD very closely and was therefore 
a more precise implementation of the desired turbulence spectrum. 


To investigate the extent to which using different random number generators affected 
the resulting motion, the vertical gusts were used in the aircraft flight dynamics simulation. 
The vertical gusts can be transformed into an angle of attack gusts as 


Aw, 
Aa, = arctan (=) (17) 


for small perturbations about the reference condition. The relative angle of attack Aa — 
Aa, manifests in the aerodynamic forces and moments [22] and the short period model is 
augmented as 


Pa ae ieee ee ee Dis ZV NS 
o)a[ae Me [at] [ a Tie] [se] 88 
Aa 1 0 0 0 
KGS) 0° | - +{ 0 0 | a (18b) 
0 Vo Vo Vo 
Aa; 9 La Go 24 a 45, ——G 4a g 


This model was simulated using each generated white noise sequence. To create variabil- 
ity in the initial conditions of the state variables, the turbulence sequences were run through 
the model with zero elevator deflections. Monte Carlo results of the vertical accelerometer 
perturbations for each different generator are shown in Fig. 13. The thick lines show the 
ensemble means for the simulations, whereas the thin lines show two standard deviations of 
scatter in the time histories. Similar trends were observed for the angle of attack and pitch 
rate outputs. The spread in the time histories were similar for each random number genera- 
tor. Although sequences from MATLAB and Random.org had less precise implementations 
of the turbulence spectrum than Fourier synthesis, as shown in Fig. 12, these errors were 
averaged out over many Monte Carlo simulations. 
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7 Concluding Remarks 


Three random number generators, used to create sequences of Gaussian white noise, were 
compared. The first generator was the MATLAB implementation of the Mersenne-Twister 
algorithm, which was included because of its widespread use. The second generator was a 
website called Random.org, which processes measurements of atmospheric noise to create 
random sequences. This generator was included because it claims to be a true random 
number generator, creating the number sequences from measurements of physical processes 
considered to be random. The third method was the method of Fourier synthesis by Lanczos 
and Bellai. This method constructs noise from definitions of amplitude and phase spectra, 
and was included because of its excellent advertised properties. 


Ensembles of 200 sample records, each having 601 numbers, were obtained for each of 
the three random number generators. These sequences were analyzed individually and as an 
ensemble, and were used in dynamic modeling applications including parameter estimation 
of stability and control derivatives and simulation of atmospheric turbulence. The same 
sample records were used throughout this paper. 


The main conclusions of this report can be summarized as follows: 


1. Each random number generator was generally of good quality and sufficient for Monte 
Carlo simulation in dynamic modeling problems. 


2. The Fourier synthesis method was consistently more accurate in approximating Gaus- 
sian white noise sequences, per the statistical tools applied. 


The deficiencies of the generators exhibited in any single sample record were generally 
averaged out over many Monte Carlo simulations, and therefore each generator is suitable 
for Monte Carlo analysis in aircraft dynamic modeling applications. All three generators 
performed similarly on average, where the differences were negligible for the applications 
considered. However, the Fourier synthesis method was significantly more consistent in 
achieving Gaussian white noise sequences than the other two methods. This was true in 
terms of the sequence mean and variance, normality, autocorrelation, and power spectral 
density, all of which are tools relevant to the analysis of aircraft dynamic modeling. 


It is therefore recommended to use Fourier synthesis for generating Gaussian white noise 
sequences in Monte Carlo simulation. These sequences required negligibly more time than 
the built-in MATLAB generator, and produced a better approximation to the white noise 
sequences assumed in the analysis. Although the examples shown here are themselves 
engineering approximations (e.g., measurement noise is not exactly white, stability and 
control derivatives are Taylor series approximations, and turbulence can take on a different 
character than the models predict), having better analysis tools that work as expected can 
simplify software development and possibly lead to fewer needed Monte Carlo simulations. 
More importantly, the analysis can be more precise and therefore the decisions based on that 
analysis will be clearer, more accurate, and more reliable, because of the greater precision 
in the white noise generation. 
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A -— Variance of the Finite Fourier Series 


For ideal white noise, selecting c, = ,/2/M results in a unit-variance sequence. This 
can be seen by considering the variance of the continuous, zero-mean variable v(t) over the 


domain ¢ € [0,T] 


2 
dt 


1 T 


where w;, = 27k/T is the radian frequency. 


M 
S- Cr Sin (wet + bx) 
k=1 


Expanding the squared expression results in terms proportional to sin?(w,z¢ + @,) and 
sin(w,t + d,) sin(wmt + dm). These latter terms sum to zero because the sinusoids are 
harmonics, which are mutually orthogonal, regardless of the phase angles. This simplifies 
the variance to 


T M 
var[u(t)] = a | S- cf sin? (wt + b,)dt 


T 
0 k=1 


which, due to linearity, can be rearranged as 


eae 
var[u(t)] = T S- a | sin? (pt + o,)dt 
k=1 


Employing double-angle identities, the above integral evaluates to 


Tr P 
1 
i sin? (wyt + dx)dt = = i 1 — cos(2wzt + 2h, )dt 
0 0 


2 

T sin(2w,T + 20x) sin(2¢z) 
= + 

2: Aw, Aw, 
= T/2 


because k is an integer and the sinusoids are periodic in T. 


Substituting this result back into the variance, 


For white noise, cz is the same for each harmonic frequency. For unit variance white noise, 
we require var[v(t)] = 1 so that 
2 


Ch = uM 


It is interesting that the orthogonality of the Fourier series means that the variance of the 
signal is only dependent on the number of frequencies included, and not upon the random 
phase angles that create constructive and destructive interference of the sinusoids. 
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B  — Autocorrelation of the Finite Fourier Series 


The autocorrelation of the Fourier synthesis sequences were good approximations of 
white noise. The reason for this can be seen by extrapolation from the sum of the first two 
harmonic sinusoids, 


u(ti) = v1 (ti) + vo(ta) 
. (2a . (An 
=c, sin (Fs + 61) + c2s1n (Fs + 6a) 


The autocorrelation function of this sequence can be expanded using the definition of 
the autocorrelation, and then rearranged as 


yy (l) =F [u(ti)] [v(t + 1)] 


= E [vi (ti) + vo(ta)] [vr (ti + 0) + ve(ts + 0] 
= Plorv1 (J) + Pougv2 (1) + Poy v2 (J) TF Pvguy (1) 


[v1 


which is the sum of the autocorrelations and crosscorrelations of the first two harmonic 
sinusoids. By extrapolation of this expression, the autocorrelation of the finite Fourier series 
is equivalent to the sum of the autocorrelations and crosscorrelations of all the sinusoids. 


The autocorrelation r,,,, and the crosscorrelation r,,,, for a single sample record ob- 
tained using Fourier synthesis is shown in Fig. 14. The autocorrelation has a coherent peak 
at lag index 1 = 0, and decaying oscillations elsewhere. The crosscorrelation is zero at 1 = 0, 
and has oscillations at other lag indices. When many sinusoids are used in the finite Fourier 
series, the peak at | = 0 is reinforced whereas the autocorrelation at the other lag indices 
averages to nearly zero. This effect creates a good approximation to the autocorrelation of 
a white noise sequence. 


21 


Tables 


Table 1: Summary of single sample records 


Characteristic Statistic Matlab Random.org Fourier synthesis 

Time history Mean +0.0189 +0.0065 +0.0009 
Variance 0.9899 1.1085 1.0000 

Normality R? with N(0, 1) 0.9960 0.9962 0.9983 
Autocorrelation Variance fy, (0) 0.9886 1.1067 0.9983 
Percentage outside 20 0.4992 2.1631 0.0000 
Power spectral density Mean 0.9926 1.1122 0.9991 

Variance 0.7887 1.1811 0.0010 


Table 2: Summary of ensembles of sample records 


Characteristic Statistic Matlab Random.org Fourier Synthesis 
Time history Average of mean —0.0014 —0.0002 +0.0002 
Scatter of mean 0.0421 0.0420 0.0015 
Average of variance 0.9969 0.9957 1.0000 
Scatter of variance 0.0544 0.0549 0.0000 
Normality Average of R? with V(0,1) 0.9948 0.9948 0.9972 
Scatter of R? with N(0, 1) 0.0029 0.0031 0.0013 
Autocorrelation Average of variance fy, (0) 0.9970 0.9958 0.9983 
Scatter of variance #,.,(0) 0.0544 0.0547 0.0000 
Percentage outside 20 0.0048 0.0032 0.0000 
Power spectral density Average of mean 0.9990 0.9972 1.0030 
Scatter of mean 0.0712 0.0679 0.0035 
Average of variance 0.9991 1.0018 0.0028 
Scatter of variance 0.1940 0.1873 0.0002 
Pairwise correlations Average of correlation —0.0006 +0.0001 +0.0001 
Scatter of correlation 0.0411 0.0408 0.0411 
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Table 3: Simulation model parameters for the T-2 aircraft 


Parameter Value Unit 
Vo 128.91 ft/s 
Zo —2.5394 1/s 
Zq —0.0551 - 
Zo. —0.2543  1/s 
Ma —36.366 1/s? 
My —3.1959 1/s 
Ms, —39.112 1/s? 
Ls 500 ft 
Ow 16.878 ft/s 


Table 4: Summary for Parameter estimation 


Parameter Statistic Matlab  Random.org Fourier synthesis 

Le Average of estimate —2.5390 —2.5403 —2.5405 
Scatter of estimate 0.0230 0.0237 0.0238 
Average of standard error 0.0212 0.0213 0.0217 
Scatter of standard error 0.0033 0.0035 0.0021 

Zq Average of estimate —0.0557 —0.0552 —0.0549 
Scatter of estimate 0.0054 0.0055 0.0046 
Average of standard error 0.0048 0.0048 0.0049 
Scatter of standard error 0.0007 0.0007 0.0004 

LZbe Average of estimate —0.2569 —0.2537 —0.2578 
Scatter of estimate 0.0254 0.0257 0.0229 
Average of standard error 0.0231 0.0231 0.0235 
Scatter of standard error 0.0036 0.0034 0.0020 

Ma Average of estimate —35.8216 —35.8363 —35.8913 
Scatter of estimate 0.3226 0.3284 0.2942 
Average of standard error 0.4047 0.4064 0.4000 
Scatter of standard error 0.0963 0.0829 0.0877 

My Average of estimate —3.0624 —3.0560 —3.0721 
Scatter of estimate 0.1150 0.1061 0.1110 
Average of standard error 0.1195 0.1192 0.1186 
Scatter of standard error 0.0235 0.0195 0.0175 

M3, Average of estimate —37.8299 —37.8128 —37.8970 
Scatter of estimate 0.5992 0.5944 0.6122 
Average of standard error 0.6225 0.6179 0.6160 
Scatter of standard error 0.1283 0.1087 0.0898 
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Figure 1: Single sample records from each random number generator 
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Figure 2: Normal probability plots, single sample records 
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Figure 3: Autocorrelations, single sample records 
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Figure 4: Power spectral densities, single sample records 
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Figure 5: Cross plots of two sample records 
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Figure 6: Efficiency of the Fourier synthesis generator in achieving a Gaussian distribution 
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Figure 7: T-2 airplane (credit: NASA Langley Research Center) 
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Figure 8: Simulated data for parameter estimation 
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Figure 9: Model fits to noisy vertical accelerometer data, single sample records 
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Figure 10: Model fits to noisy pitch acceleration data, single sample records 
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Figure 11: Vertical gust time histories, single sample records 
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Figure 12: Vertical gust power spectral densities, single sample records 
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Figure 13: Simulated vertical accelerations with turbulence, Monte Carlo results 
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Figure 14: Autocorrelation and crosscorrelation of harmonic sinusoids 
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